# Seurat
library(Seurat)
# Single R
library(SingleCellExperiment)
library(SingleR)
library(celldex)
# Data
library(dplyr)
library(Matrix)
# Plotting
library(ggplot2)
library(RColorBrewer)
library(patchwork)
Attaching SeuratObject
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
colWeightedMeans, colWeightedMedians, colWeightedSds,
colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
rowWeightedSds, rowWeightedVars
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: ‘BiocGenerics’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, append, as.data.frame, basename, cbind, colnames,
dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unsplit, which.max, which.min
Loading required package: S4Vectors
Attaching package: ‘S4Vectors’
The following objects are masked from ‘package:base’:
expand.grid, I, unname
Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
Attaching package: ‘SummarizedExperiment’
The following object is masked from ‘package:SeuratObject’:
Assays
The following object is masked from ‘package:Seurat’:
Assays
Attaching package: ‘celldex’
The following objects are masked from ‘package:SingleR’:
BlueprintEncodeData, DatabaseImmuneCellExpressionData,
HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData,
MouseRNAseqData, NovershternHematopoieticData
Attaching package: ‘dplyr’
The following object is masked from ‘package:Biobase’:
combine
The following objects are masked from ‘package:GenomicRanges’:
intersect, setdiff, union
The following object is masked from ‘package:GenomeInfoDb’:
intersect
The following objects are masked from ‘package:IRanges’:
collapse, desc, intersect, setdiff, slice, union
The following objects are masked from ‘package:S4Vectors’:
first, intersect, rename, setdiff, setequal, union
The following objects are masked from ‘package:BiocGenerics’:
combine, intersect, setdiff, union
The following object is masked from ‘package:matrixStats’:
count
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Attaching package: ‘Matrix’
The following object is masked from ‘package:S4Vectors’:
expand
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
# Filtering Parameter
nFeature_RNA_min_m <- 500
nFeature_RNA_min_p <- 500
nFeature_RNA_max_m <- 3000
nFeature_RNA_max_p <- 3000
nCount_RNA_min_m <- 1500
nCount_RNA_min_p <- 1500
nCount_RNA_max_m <- 20000
nCount_RNA_max_p <- 20000
pMt_RNA_max_m <- 7.5
pMt_RNA_max_p <- 7.5
# Files
so_raw_file <- "data/object/raw.rds"
so_qc_file <- "data/object/seurat.rds"
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so_raw <- readRDS(so_raw_file)
so_raw$nFeature_RNA_max <- ifelse(so_raw$tissue == "Myeloid", nFeature_RNA_max_m, nFeature_RNA_max_p)
so_raw$nFeature_RNA_min <- ifelse(so_raw$tissue == "Myeloid", nFeature_RNA_min_m, nFeature_RNA_min_p)
so_raw$nCount_RNA_max <- ifelse(so_raw$tissue == "Myeloid", nCount_RNA_max_m, nCount_RNA_max_p)
so_raw$nCount_RNA_min <- ifelse(so_raw$tissue == "Myeloid", nCount_RNA_min_m, nCount_RNA_min_p)
so_raw$pMt_RNA_max <- ifelse(so_raw$tissue == "Myeloid", pMt_RNA_max_m, pMt_RNA_max_p)
# QC
so_raw$qc_class <- ifelse(
so_raw$cellranger_class == "Cell" &
so_raw$nFeature_RNA <= so_raw$nFeature_RNA_max &
so_raw$nFeature_RNA >= so_raw$nFeature_RNA_min &
so_raw$nCount_RNA <= so_raw$nCount_RNA_max &
so_raw$nCount_RNA >= so_raw$nCount_RNA_min &
so_raw$pMt_RNA <= so_raw$pMt_RNA_max,
"pass", "fail"
)
Empty droplets were determined with CellRanger V3.0.2 Lun et al., 2019 EmptyDrop heuristic. RNAse activity of granulocytes might be wrongly identified as empty cells by CellRanger.
Typical Sample A steep drop-off is indicative of good separation between the cell-associated barcodes and the barcodes associated with empty GEMs. A ideal barcode rank plot has a distincitve shape, which is referred to as a "cliff and knee".
Heterogeneous Sample Heterogeneous populations of cells in a sample result in two "cliff and knee" distributions. However, there should still be clear separation between the bacodes.
Compromised Sample Round curve and lack of steep cliff may indicate low sample quality or loss of single-cell behavior. This can be due to a wetting failure, premature cell lysis, or low cell viability.
Compromised Sample Defined cliff and knee, but the total number of barcodes detected may be lower than expected. This can be caused by a sample clog or inaccurate cell count.
rank_plot <- ggplot(so_raw@meta.data, aes(x = log10(nCount_RNA_rank), y = log10(nCount_RNA), color = cellranger_class)) +
geom_point() +
geom_hline(aes(yintercept = log10(nCount_RNA_min)), color = "red", linetype = "longdash") +
scale_color_manual(values = color$cellranger_class) +
ggtitle("Barcode rank plot") +
xlab("log10(cell barcode rank)") + ylab("log10(cell UMI counts)") +
facet_grid(tissue~treatment) +
theme(aspect.ratio = 1, legend.position = "bottom")
options(repr.plot.width = 5, repr.plot.height = 5)
rank_plot
ggsave(rank_plot, filename = "result/plot/seurat/rank_plot.png", width = 4, height = 4)
so_qc <- subset(so_raw, subset = cellranger_class == "Cell")
so_qc$sample_rep <- gsub("^.*_", "", so_qc$sample_name)
qc_1 <- ggplot(so_qc@meta.data, aes(x = log10(nCount_RNA), fill = tissue)) +
geom_density() +
ggtitle("Density plot UMI count") + xlab("log10(UMI count)") + ylab("Density") +
geom_vline(aes(xintercept = log10(nCount_RNA_min)), color = "red", linetype = "longdash") +
geom_vline(aes(xintercept = log10(nCount_RNA_max)), color = "red", linetype = "longdash") +
scale_x_continuous(breaks = integer_breaks()) +
scale_fill_manual(values = color$tissue) +
facet_grid(tissue~treatment+sample_rep) +
theme(legend.position = "bottom", aspect.ratio = 1)
qc_2 <- ggplot(so_qc@meta.data, aes(x = log10(nFeature_RNA), fill = tissue)) +
geom_density() +
ggtitle("Density plot Feature count") + xlab("log10(Feature count)") + ylab("Density") +
geom_vline(aes(xintercept = log10(nFeature_RNA_min)), color = "red", linetype = "longdash") +
geom_vline(aes(xintercept = log10(nFeature_RNA_max)), color = "red", linetype = "longdash") +
scale_x_continuous(breaks = integer_breaks()) +
scale_fill_manual(values = color$tissue) +
facet_grid(tissue~treatment+sample_rep) +
theme(legend.position = "bottom", aspect.ratio = 1)
qc_3 <- ggplot(so_qc@meta.data, aes(x = pMt_RNA, fill = tissue)) +
geom_density() +
ggtitle("Density plot Mt %") + xlab("Mt [%]") + ylab("Density") +
geom_vline(aes(xintercept = pMt_RNA_max), color = "red", linetype = "longdash") +
scale_x_continuous(breaks = integer_breaks()) +
scale_fill_manual(values = color$tissue) +
facet_grid(tissue~treatment+sample_rep) +
theme(legend.position = "bottom", aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 10)
qc_1 + qc_2 + qc_3 + plot_layout(ncol = 2) & theme(legend.position = "bottom")
ggsave(qc_1, filename = "result/plot/seurat/density_umi.png", width = 4, height = 4)
ggsave(qc_2, filename = "result/plot/seurat/density_feature.png", width = 4, height = 4)
ggsave(qc_3, filename = "result/plot/seurat/density_mt.png", width = 4, height = 4)
sc_1 <- ggplot(so_qc@meta.data, aes(x = log10(nCount_RNA), y = log10(nFeature_RNA), color = pMt_RNA)) +
geom_point() + ggtitle("Mitochondrial gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") +
geom_vline(aes(xintercept = log10(nCount_RNA_min)), color = "red", linetype = "longdash") +
geom_vline(aes(xintercept = log10(nCount_RNA_max)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_min)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_max)), color = "red", linetype = "longdash") +
facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio = 1, legend.position = "bottom") +
scale_size(guide = guide_legend(direction = "vertical"))
sc_2 <- ggplot(so_qc@meta.data, aes(x = log10(nCount_RNA), y = log10(nFeature_RNA), color = pHb_RNA)) +
geom_point() + ggtitle("Hemoglobin gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") +
geom_vline(aes(xintercept = log10(nCount_RNA_min)), color = "red", linetype = "longdash") +
geom_vline(aes(xintercept = log10(nCount_RNA_max)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_min)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_max)), color = "red", linetype = "longdash") +
facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio = 1, legend.position = "bottom") +
scale_size(guide = guide_legend(direction = "vertical"))
sc_3 <- ggplot(so_qc@meta.data, aes(x = log10(nCount_RNA), y = log10(nFeature_RNA), color = pRp_RNA)) +
geom_point() + ggtitle("Ribsonmal gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") +
geom_vline(aes(xintercept = log10(nCount_RNA_min)), color = "red", linetype = "longdash") +
geom_vline(aes(xintercept = log10(nCount_RNA_max)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_min)), color = "red", linetype = "longdash") +
geom_hline(aes(yintercept = log10(nFeature_RNA_max)), color = "red", linetype = "longdash") +
facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio = 1, legend.position = "bottom") +
scale_size(guide = guide_legend(direction = "vertical"))
options(repr.plot.width = 15, repr.plot.height = 10)
sc_1 + sc_2 + sc_3 + plot_layout(ncol = 2) & theme(legend.position = "bottom")
ggsave(sc_1, filename = "result/plot/seurat/sc_mt.png", width = 4, height = 4)
ggsave(sc_2, filename = "result/plot/seurat/sc_hg.png", width = 4, height = 4)
ggsave(sc_3, filename = "result/plot/seurat/sc_rb.png", width = 4, height = 4)
so_qc <- subset(so_qc, subset = nFeature_RNA >= nFeature_RNA_min & nFeature_RNA <= nFeature_RNA_max & nCount_RNA >= nCount_RNA_min & nCount_RNA <= nCount_RNA_max & pMt_RNA <= pMt_RNA_max)
cc_kowalczyk <- read.csv("cc_kowalczyk.csv")
cc_kowalczyk <- cc_kowalczyk[cc_kowalczyk$sig_population >= 5, ]
cc_kowalczyk <- cc_kowalczyk[cc_kowalczyk$gene %in% rownames(so_raw), ]
cc_kowalczyk_g0 <- cc_kowalczyk[cc_kowalczyk$state == "down", ]
cc_kowalczyk_int <- cc_kowalczyk[cc_kowalczyk$state == "up", ]
cc_kowalczyk_s <- cc_kowalczyk_int[cc_kowalczyk_int$S > cc_kowalczyk_int$G2_M, ]
cc_kowalczyk_g2m <- cc_kowalczyk_int[cc_kowalczyk_int$G2_M > cc_kowalczyk_int$S, ]
so_raw <- CellCycleScoring(so_raw, s.features = cc_kowalczyk_s$gene, g2m.features = cc_kowalczyk_g2m$gene, set.ident = FALSE)
colnames(so_raw@meta.data) <- gsub("Phase", "cc_phase_class", colnames(so_raw@meta.data))
colnames(so_raw@meta.data) <- gsub("S.Score", "msS_RNA", colnames(so_raw@meta.data))
colnames(so_raw@meta.data) <- gsub("G2M.Score", "msInt_RNA", colnames(so_raw@meta.data))
so_raw$msCC_diff_RNA <- so_raw$msS_RNA - so_raw$msInt_RNA
so_raw <- AddModuleScore(so_raw, features = list(cc_kowalczyk_g0$gene), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msG0_RNA")
so_raw <- AddModuleScore(so_raw, features = list(cc_kowalczyk_int$gene), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msInt_RNA")
qc_vln_FUN <- function(data, y, fill, ylab = "", ymin, ymax) {
vln_plot_1 <- ggplot(data, aes(x = treatment, y = {{y}}, color = {{fill}})) +
geom_jitter(alpha = 0.2, shape = 16, color = "gray") +
geom_boxplot(alpha = 1.0) + xlab("") +
ylim(ymin, ymax) +
scale_color_manual(values = color$tissue) +
ggtitle("Cell containing GEM") +
facet_wrap(~tissue, scales = "free_x") + ylab(ylab) +
theme(
plot.title = element_text(size = 12, face = "bold", margin = margin(t = 0, r = 0, b = 5, l = 0)),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
strip.text = element_blank()
)
vln_plot_2 <- ggplot(data[data$qc_class == "pass", ], aes(x = treatment, y = {{y}}, color = {{fill}})) +
geom_jitter(alpha = 0.2, shape = 16, color = "gray") +
geom_boxplot(alpha = 1.0) + xlab("") +
ylim(ymin, ymax) +
scale_color_manual(values = color$tissue) +
ggtitle("Filtered") +
facet_wrap(~tissue, scales = "free_x") + ylab(ylab) +
theme(
plot.title = element_text(size = 12, face = "bold", margin = margin(t = 0, r = 0, b = 5, l = 0)),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
strip.text = element_blank()
)
vln_plot <- vln_plot_1 + vln_plot_2 + plot_layout(guides = "collect") & theme(legend.position = "bottom")
return(vln_plot)
}
data <- dplyr::filter(so_raw@meta.data, cellranger_class == "Cell")
qc_vln_1 <- qc_vln_FUN(data, nCount_RNA, tissue, "UMI [count]", ymin = 0, ymax = max(so_raw$nCount_RNA)+100)
qc_vln_2 <- qc_vln_FUN(data, nFeature_RNA, tissue, "Feature [count]", ymin = 0, ymax = max(so_raw$nFeature_RNA)+100)
qc_vln_3 <- qc_vln_FUN(data, msCC_diff_RNA, tissue, "cc diff [count]", ymin = min(so_raw$msCC_diff_RNA), ymax = max(so_raw$msCC_diff_RNA))
qc_vln_4 <- qc_vln_FUN(data, pMt_RNA, tissue, "Mt [%]", ymin = 0, ymax = 100)
qc_vln_5 <- qc_vln_FUN(data, pHb_RNA, tissue, "Hb [%]", ymin = 0, ymax = 100)
qc_vln_6 <- qc_vln_FUN(data, pRp_RNA, tissue, "Rbl [%]", ymin = 0, ymax = 100)
options(repr.plot.width = 10, repr.plot.height = 5)
qc_vln_1
qc_vln_2
qc_vln_3
qc_vln_4
qc_vln_5
qc_vln_6
ggsave(qc_vln_1, filename = "result/plot/seurat/qc_vln_1.png", width = 5, height = 2.5)
ggsave(qc_vln_2, filename = "result/plot/seurat/qc_vln_2.png", width = 5, height = 2.5)
ggsave(qc_vln_3, filename = "result/plot/seurat/qc_vln_3.png", width = 5, height = 2.5)
ggsave(qc_vln_4, filename = "result/plot/seurat/qc_vln_4.png", width = 5, height = 2.5)
ggsave(qc_vln_5, filename = "result/plot/seurat/qc_vln_5.png", width = 5, height = 2.5)
ggsave(qc_vln_6, filename = "result/plot/seurat/qc_vln_6.png", width = 5, height = 2.5)
Warning message: “Removed 1 rows containing missing values (geom_point).” Warning message: “Removed 1 rows containing missing values (geom_point).”
Warning message: “Removed 29 rows containing missing values (geom_point).” Warning message: “Removed 1 rows containing missing values (geom_point).”
Warning message: “Removed 6192 rows containing missing values (geom_point).” Warning message: “Removed 3026 rows containing missing values (geom_point).”
Warning message: “Removed 3 rows containing missing values (geom_point).”
Warning message: “Removed 22 rows containing missing values (geom_point).” Warning message: “Removed 1 rows containing missing values (geom_point).” Warning message: “Removed 6274 rows containing missing values (geom_point).” Warning message: “Removed 3014 rows containing missing values (geom_point).” Warning message: “Removed 6 rows containing missing values (geom_point).”
so_qc <- CellCycleScoring(so_qc, s.features = cc_kowalczyk_s$gene, g2m.features = cc_kowalczyk_g2m$gene, set.ident = FALSE)
colnames(so_qc@meta.data) <- gsub("Phase", "cc_phase_class", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("S.Score", "msS_RNA", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("G2M.Score", "msInt_RNA", colnames(so_qc@meta.data))
so_qc$msCC_diff_RNA <- so_qc$msS_RNA - so_qc$msInt_RNA
so_qc <- AddModuleScore(so_qc, features = list(cc_kowalczyk_g0$gene), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msG0_RNA")
so_qc <- AddModuleScore(so_qc, features = list(cc_kowalczyk_int$gene), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msInt_RNA")
so_qc <- NormalizeData(
object = so_qc,
assay = "RNA",
normalization.method = "LogNormalize",
scale.factor = 10000
)
so_qc <- FindVariableFeatures(
object = so_qc,
assay = "RNA",
selection.method = "vst",
nfeatures = 3000
)
variable_features <- VariableFeatures(so_qc)
VariableFeatures(so_qc) <- variable_features[!variable_features %in% c(cc_kowalczyk$gene)]
Regress nCount_RNA (total UMI per cell) and pMT_RNA (mitochondrial UMI percentage)
so_qc <- ScaleData(
object = so_qc,
assay = "RNA",
vars.to.regress = c("nCount_RNA", "pMt_RNA", "msCC_diff_RNA"),
do.scale = TRUE,
do.center = TRUE
)
Regressing out nCount_RNA, pMt_RNA, msCC_diff_RNA Centering and scaling data matrix
source("bin/seurat_function.R")
so_qc <- dim_clust(
# Seurat
so = so_qc,
assay = "RNA",
# PCA features
blacklist_genes = NULL,
# Dim reduction
dims_pca = 100, # Default 20 - RunPCA dims
dims_umap = 90, # Default NULL - RunUMAP dims
dims_tsne = 30, # Default 1:5
min_dist = 0.3, # Default 0.3 - RunUMAP dmin.ist - controls how tightly the embedding
# Cluster
dims_cluster = 10, # Default 1:10 - FindNeighbors dims
cluster_res = 0.8 # Default 0.8 - FindClusters resoluton - above(below) 1.0 to obtain larger(smaller) number of communities)
)
PC_ 1 Positive: Car2, Prdx2, Blvrb, Rps2, Glrx5, Hbb-bs, Hba-a1, Hbb-bt, Cd24a, Car1 Ctse, Cpox, Hba-a2, Rhd, Ybx3, Hsp90aa1, Mgst3, Eif5a, Odc1, C1qtnf12 Hdgf, Ermap, Hemgn, Cldn13, Gypa, Ncl, Stard10, Klf1, Rps17, Ubac1 Negative: Tyrobp, Tmsb4x, Ptprc, Laptm5, Cd52, Fcer1g, Ctss, Cd74, Cyba, Coro1a Lcp1, Lsp1, Cst3, Sh3bgrl3, Napsa, Gm2a, Rgs2, Lst1, Psap, Sat1 Spi1, Arpc1b, Gngt2, H2-Aa, Fyb, Pld4, Pou2f2, Cybb, H2-Eb1, H2-Ab1 PC_ 2 Positive: Ciita, Dnase1l3, H2-Eb1, H2-Ab1, Ccnd1, H2-Aa, H2-Oa, S100a11, Slamf7, Wdfy4 Tbc1d4, Xcr1, H2-DMa, Jaml, Batf3, H2-DMb2, Ptms, Flt3, Ppp1r14a, Adam23 Ffar2, Cd8a, Cd74, Avpi1, Siglecg, H2-Ob, Rab7b, Clec9a, H2-DMb1, Gsn Negative: Csf1r, Cebpb, Adgre4, Ace, Ctsb, Lrp1, Fcgr4, Pilra, Clec4a3, Lyz2 Msrb1, Hip1, Cd300e, Cd300ld, Cx3cr1, Fabp4, Hpgd, Smpdl3a, Cyp4f18, Pou2f2 Ceacam1, Itgal, Tnfrsf1b, Pilrb2, Stk10, Treml4, Cybb, Nupr1, Tppp3, Smpdl3b PC_ 3 Positive: Ace, Itgal, Pglyrp1, Spn, Lyz2, Cyp4f18, Eno3, Plac8, Cx3cr1, Msrb1 Ceacam1, Ifitm2, Klf4, Klf2, Nr4a1, Emp3, Cebpb, Ifitm6, Tppp3, Coro1a Smpdl3b, Hp, Nupr1, Il17ra, Stap1, Rbpms, Cd300ld, Ifitm3, Ear2, Myo1g Negative: Vcam1, C1qc, C1qa, C1qb, Fcna, Cd5l, Mertk, Igf1, Mrc1, Axl Maf, Itgad, Slc40a1, Mafb, Spic, Ccr3, Hmox1, Frmd4b, Ptgs1, Epb41l3 Slpi, Kcna2, Kcnj10, Hpgds, Sema6d, Gfra2, Sdc3, Dst, Pld3, Clec4n PC_ 4 Positive: Fth1, Hbb-bt, Hba-a1, Hba-a2, Hbb-bs, Ctse, Cd24a, Gypa, Mgst3, Rhd Cldn13, Ftl1, Slc4a1, Hebp1, Hemgn, Cpox, Tspan33, Tfrc, Ubac1, Hagh Smim1, Prdx2, Isg20, Glrx5, Blvrb, Slc25a37, Car2, Carhsp1, Rhag, Lmna Negative: Cmtm7, S100a10, Myb, Ptprcap, Gata2, Tagln2, Sox4, Cd63, Gm15915, Cpa3 Lgals1, Apoe, Ctsg, Ifitm1, Ctla2a, Ms4a2, Vim, Serpinb1a, Srgn, Phgdh Angpt1, Ifitm2, Mcpt8, Muc13, Ikzf2, Fcer1a, Gmfg, Plac8, Ms4a3, Cst7 PC_ 5 Positive: Ccl5, S100a4, Ccnd1, Relb, Ffar2, Siglecg, Dtx1, Ltb, Tbc1d4, Tmem176b Cd7, H2-DMb2, Il4i1, Plxnc1, Tmem176a, Ppp1r14a, Cyp4f16, Rgs1, Traf1, Lyst Ryr3, H2-Eb2, Arap2, St8sia6, Avpi1, Tctex1d2, Man1a, Icosl, Ube2e2, Adam23 Negative: Xcr1, Naaa, Rab7b, Irf8, Cd8a, Clec9a, Ppt1, Tlr3, Cadm1, A530099J19Rik P2ry14, Mpeg1, AC163354.1, Ifi205, Pianp, Itgae, Pdia5, Gcsam, Lgals3, Olfm1 Sept3, Cst3, Aif1, Tlr11, Dbn1, Crip1, Btla, Cxcl9, Slc8b1, Ece1 Computing nearest neighbor graph Computing SNN Computing nearest neighbors Only one graph name supplied, storing nearest-neighbor graph only Warning message: “The following arguments are not used: assay” Warning message: “The following arguments are not used: assay”
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 22201 Number of edges: 695360 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.8892 Number of communities: 22 Elapsed time: 3 seconds
Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 14:15:25 UMAP embedding parameters a = 0.9922 b = 1.112 14:15:26 Commencing smooth kNN distance calibration using 1 thread 14:15:29 Initializing from normalized Laplacian + noise 14:15:31 Commencing optimization for 200 epochs, with 565500 positive edges 14:15:54 Optimization finished 14:15:54 UMAP embedding parameters a = 0.9922 b = 1.112 14:15:54 Read 22201 rows and found 90 numeric columns 14:15:54 Using Annoy for neighbor search, n_neighbors = 20 14:15:54 Building Annoy index with metric = euclidean, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 14:15:57 Writing NN index file to temp file /tmp/Rtmp9uQVAH/file226c1df1d2d6 14:15:57 Searching Annoy index using 1 thread, search_k = 2000 14:16:02 Annoy recall = 100% 14:16:04 Commencing smooth kNN distance calibration using 1 thread 14:16:07 Initializing from normalized Laplacian + noise 14:16:08 Commencing optimization for 200 epochs, with 687332 positive edges 14:16:34 Optimization finished Warning message: “Cannot add objects with duplicate keys (offending key: UMAP_), setting key to 'rna_umap_'”
SingleR identifies marker genes from the reference and uses them to compute assignment scores (based on the Spearman correlation across markers) for each cell in the test dataset against each label in the reference. The label with the highest score is the assigned to the test cell, possibly with further fine-tuning to resolve closely related labels.
first.labels: Labels before fine-tuning
labels: Labels after fine-tuning
pruning: labels after pruning
ref <- readRDS("data/haemosphere/se_haemosphere.rds")
# Seurat to SingleCellExperiment
sce <- SingleCellExperiment(list(counts = so_qc@assays$RNA@counts))
# Predict labels
label_main <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.main, assay.type.test = "counts", de.method = "classic")
label_fine <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.fine, assay.type.test = "counts", de.method = "classic")
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(main_labels = pruned.labels, main_delta_score = tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(fine_labels = pruned.labels, fine_delta_score = tuning.scores.first)
so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))
# Set factor level for labels
so_qc$main_labels <- factor(so_qc$main_labels, levels = names(color$main_labels_haemosphere))
so_qc$fine_labels <- factor(so_qc$fine_labels, levels = names(color$fine_labels_haemosphere))
reduction <- "rna_umap_nno"
dplot_1 <- DimPlot(so_qc, reduction = reduction, group.by = "rna_snn_res.0.8", label = TRUE) &
theme(aspect.ratio = 1, legend.position = "none")
dplot_2 <- DimPlot(so_qc, reduction = reduction, group.by = "main_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so_qc, reduction = reduction, group.by = "fine_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so_qc, reduction = reduction, group.by = "tissue", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$tissue, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so_qc, reduction = reduction, group.by = "treatment", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$treatment, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
options(repr.plot.width = 15, repr.plot.height = 10)
dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + plot_layout(ncol = 3)
dplot
ggsave(dplot, filename = "result/plot/seurat/dimplot_1.png", width = 15, height = 10)
fplot_1 <- FeaturePlot(so_qc, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_qc, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_qc, reduction = reduction, features = "main_delta_score") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 5)
fplot <- fplot_1 + fplot_2 + fplot_3 + plot_layout(ncol = 3)
fplot
ggsave(fplot, filename = "result/plot/seurat/fplot_1.png", width = 9, height = 6)
fplot_1 <- FeaturePlot(so_qc, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_qc, reduction = reduction, features = "pRp_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_qc, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 5)
fplot <- fplot_1 + fplot_2 + fplot_3 + plot_layout(ncol = 3)
fplot
ggsave(fplot, filename = "result/plot/seurat/fplot_2.png", width = 9, height = 6)
fplot_1 <- FeaturePlot(so_qc, reduction = reduction, features = "msG0_RNA1") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_qc, reduction = reduction, features = "msInt_RNA1") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_qc, reduction = reduction, features = "msCC_diff_RNA") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 5)
fplot <- fplot_1 + fplot_2 + fplot_3 + plot_layout(ncol = 3)
fplot
ggsave(fplot, filename = "result/plot/seurat/fplot_3.png", width = 9, height = 6)
cluster_tissue <- ggplot(so_qc@meta.data, aes(x = seurat_clusters, fill = tissue)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$tissue) +
ggtitle("Cluster frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
cluster_treatment <- ggplot(so_qc@meta.data, aes(x = seurat_clusters, fill = treatment)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$treatment) +
ggtitle("Cluster frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
options(repr.plot.width = 10, repr.plot.height = 5)
cluster_tissue + cluster_treatment + plot_layout(ncol = 2)
ggsave(cluster_tissue, filename = "result/plot/seurat/cluster_tissue.png", width = 6, height = 3)
ggsave(cluster_treatment, filename = "result/plot/seurat/cluster_treatment.png", width = 6, height = 3)
bar_1 <- ggplot(so_qc@meta.data, aes(x = seurat_clusters, y = pHb_RNA)) + geom_boxplot()
bar_2 <- ggplot(so_qc@meta.data, aes(x = seurat_clusters, y = pMt_RNA)) + geom_boxplot()
bar_3 <- ggplot(so_qc@meta.data, aes(x = seurat_clusters, y = pRp_RNA)) + geom_boxplot()
options(repr.plot.width = 10, repr.plot.height = 10)
bar_1 + bar_2 + bar_3 + plot_layout(ncol = 2)
fplot_1 <- FeaturePlot(so_qc, reduction = reduction, features = "Ptprc") & theme(aspect.ratio = 1) & ggtitle("Ptprc (CD45)")
fplot_2 <- FeaturePlot(so_qc, reduction = reduction, features = "Cd34") & theme(aspect.ratio = 1) & ggtitle("Cd34")
fplot_3 <- FeaturePlot(so_qc, reduction = reduction, features = "Kit") & theme(aspect.ratio = 1) & ggtitle("Kit")
fplot_4 <- FeaturePlot(so_qc, reduction = reduction, features = "Tfrc") & theme(aspect.ratio = 1) & ggtitle("Tfrc (CD71)")
fplot_5 <- FeaturePlot(so_qc, reduction = reduction, features = "Gata2") & theme(aspect.ratio = 1) & ggtitle("Gata2")
fplot_6 <- FeaturePlot(so_qc, reduction = reduction, features = "Gata1") & theme(aspect.ratio = 1) & ggtitle("Gata1")
fplot <- fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot
ggsave(fplot, filename = "result/plot/seurat/fplot_2.png", width = 9, height = 6)
H2-K1 - Class I histocompatibility antigen, kappa-B alpha chain
H2-D1 - Class I histocompatibility antigen, D-B alpha chain
H2-L1 - not found
H2-Aa - Class II histocompatibility antigen, A-B alpha chain
H2-Ab1 - Class II histocompatibility antigen, A beta chain
H2-Eb1 - Class II histocompatibility antigen, I-E beta chain
H2-Eb2 - Class II histocompatibility antigen, E beta chain
mhcI_genes <- c("H2-K1", "H2-D1")
mhcII_genes <- c("H2-Aa", "H2-Ab1", "H2-Eb1", "H2-Eb2")
cd4_genes <- c("Cd4")
cd8_genes <- c("Cd8a")
so_qc <- AddModuleScore(so_qc, features = list(mhcI_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msMHCI_RNA")
so_qc <- AddModuleScore(so_qc, features = list(mhcII_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msMHCII_RNA")
so_qc <- AddModuleScore(so_qc, features = list(cd4_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msCd4_RNA")
so_qc <- AddModuleScore(so_qc, features = list(cd8_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msCd8_RNA")
fplot_mhcI <- FeaturePlot(so_qc, reduction = reduction, features = "msMHCI_RNA1") & theme(aspect.ratio = 1) & ggtitle("MHC class I")
fplot_mhcII <- FeaturePlot(so_qc, reduction = reduction, features = "msMHCII_RNA1") & theme(aspect.ratio = 1) & ggtitle("MHC class II")
fplot_cd4 <- FeaturePlot(so_qc, reduction = reduction, features = "msCd4_RNA1") & theme(aspect.ratio = 1) & ggtitle("CD4")
fplot_cd8 <- FeaturePlot(so_qc, reduction = reduction, features = "msCd8_RNA1") & theme(aspect.ratio = 1) & ggtitle("CD8")
fplot_mhc <- fplot_mhcI + fplot_mhcII + fplot_cd4 + fplot_cd8 + plot_layout(ncol = 2)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot_mhc
ggsave(fplot_mhc, filename = "result/plot/seurat/fplot_mhc.png", width = 6, height = 6)
Trac - T cell receptor alpha constant
Trbc1 - T cell receptor beta constant 1
Trbc2 - T cell receptor beta constant 2
Trdc - T cell receptor delta constant
Trgc1 - T cell receptor gamma constant 1
Trgc2 - T cell receptor gamma constant 2
Cd247 - T-cell surface glycoprotein Cd3 zeta chain (Cd3z)
Cd3g - T-cell surface glycoprotein Cd3 gamma chain
Cd3e - T-cell surface glycoprotein Cd3 epsilon chain
Cd3d - T-cell surface glycoprotein Cd3 delta chain
tcr_genes <- c("Trac", "Trbc1", "Trbc2", "Trdc", "Trgc1", "Trgc2")
tcr_cd3_genes <- c("Cd247", "Cd3g", "Cd3e", "Cd3d")
so_qc <- AddModuleScore(so_qc, features = list(tcr_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msTcr_RNA")
so_qc <- AddModuleScore(so_qc, features = list(tcr_cd3_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msTcr_cd3_RNA")
fplot_tcr <- FeaturePlot(so_qc, reduction = reduction, features = "msTcr_RNA1") & theme(aspect.ratio = 1) & ggtitle("TCR")
fplot_tcr_cd3 <- FeaturePlot(so_qc, reduction = reduction, features = "msTcr_cd3_RNA1") & theme(aspect.ratio = 1) & ggtitle("Cd247/Cd3 family")
fplot_tcr_cd3 <- fplot_tcr + fplot_tcr_cd3 + plot_layout(ncol = 2)
options(repr.plot.width = 10, repr.plot.height = 5)
fplot_tcr_cd3
ggsave(fplot_tcr_cd3, filename = "result/plot/seurat/fplot_tcr_cd3.png", width = 6, height = 3)
Warning message: “The following features are not present in the object: Trgc1, Trgc2, not searching for symbol synonyms”
Naive B-cells produce the following Ig classes:
Ighm - Immunoglobulin heavy constant mu (naive B-cells)
Ighd - Immunoglobulin heavy constant delta (naive B-cells)
Through isotope switching the following Ig classes can be produced:
Ighg1 - Immunoglobulin heavy constant gamma (Mouse with Igh1b have Igg2c isotope instead of Igg2a)
Ighg2a - Immunoglobulin heavy constant gamma (NA)
Ighg2b - Immunoglobulin heavy constant gamma
Ighg2c - Immunoglobulin heavy constant gamma
Ighg3 - Immunoglobulin heavy constant gamma
Igha - Immunoglobulin heavy constant alpha
Ighe - Immunoglobulin heavy constant epsilon (NA)
Igkc - Immunoglobulin kappa constant (light chain)
Iglc - Immunoglobulin lambda constant (light chain)
Ighm_genes <- c("Ighm")
Ighd_genes <- c("Ighd")
Ighg_genes <- c("Ighg1", "Ighg2a", "Ighg2b", "Ighg2c", "Ighg3")
Igha_genes <- c("Igha")
Ighe_genes <- c("Ighe")
Igkc_genes <- c("Igkc")
Iglc_genes <- c("Iglc1", "Iglc2", "Iglc3")
so_qc <- AddModuleScore(so_qc, features = list(Ighm_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIghm_RNA")
so_qc <- AddModuleScore(so_qc, features = list(Ighd_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIghd_RNA")
so_qc <- AddModuleScore(so_qc, features = list(Ighg_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIghg_RNA")
so_qc <- AddModuleScore(so_qc, features = list(Igha_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIgha_RNA")
#so_qc <- AddModuleScore(so_qc, features = list(Ighe_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIghe_RNA")
so_qc <- AddModuleScore(so_qc, features = list(Igkc_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIgkc_RNA")
so_qc <- AddModuleScore(so_qc, features = list(Iglc_genes), assays = "RNA", slot = "data", ctrl = 100, nbin = 25, name = "msIglc_RNA")
fplot_Ighm <- FeaturePlot(so_qc, reduction = reduction, features = "msIghm_RNA1") & theme(aspect.ratio = 1) & ggtitle("Ighm")
fplot_Ighd <- FeaturePlot(so_qc, reduction = reduction, features = "msIghd_RNA1") & theme(aspect.ratio = 1) & ggtitle("Ighd")
fplot_Ighg <- FeaturePlot(so_qc, reduction = reduction, features = "msIghg_RNA1") & theme(aspect.ratio = 1) & ggtitle("Ighg")
fplot_Igha <- FeaturePlot(so_qc, reduction = reduction, features = "msIgha_RNA1") & theme(aspect.ratio = 1) & ggtitle("Igha")
fplot_Igkc <- FeaturePlot(so_qc, reduction = reduction, features = "msIgkc_RNA1") & theme(aspect.ratio = 1) & ggtitle("Igkc")
fplot_Iglc <- FeaturePlot(so_qc, reduction = reduction, features = "msIglc_RNA1") & theme(aspect.ratio = 1) & ggtitle("Iglc")
fplot_ig <- fplot_Ighm + fplot_Ighd + fplot_Ighg + fplot_Igha + fplot_Igkc + fplot_Iglc + plot_layout(ncol = 3)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot_ig
ggsave(fplot_ig, filename = "result/plot/seurat/fplot_ig.png", width = 9, height = 6)
Warning message: “The following features are not present in the object: Ighg2a, not searching for symbol synonyms”
# Write tsv file for each sample into 10x sample dir
for (sample_path in unique(so_qc$sample_path)) {
tsv_file <- paste0(sample_path, "/filtered_feature_bc_matrix/barcodes.qc.tsv")
print(tsv_file)
cell_id <- so_qc@meta.data[so_qc@meta.data$sample_path == sample_path, ]$cell_id
library(stringr)
write.table(x = cell_id,
file = tsv_file,
row.names = FALSE,
col.names = FALSE,
sep = "\t",
quote = FALSE)
}
[1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M1_CpG_Mac_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M1_CpG_Prog_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M2_CpG_Mac_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M2_CpG_Prog_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M6_control_Mac_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M6_control_Prog_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M8_control_Mac_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv" [1] "data/BSA_0355_SM01_10x_SPLENO/OUT/COUNT/M8_control_Prog_transcriptome/filtered_feature_bc_matrix/barcodes.qc.tsv"
# Write Seurat object
saveRDS(so_qc, so_qc_file)
Its only possible to store two out of three count matrix. Data can be recovered from counts but not vice versa. However, data can not be removed and is prioretized before counts. Therefore we fill counts in the data slot so that it gets exported.
# # Meta data
# write.csv(so_qc@meta.data, "data/object/components/meta_data/seurat_meta.csv")
# # Cluster
# write.table(so_qc[["rna_snn_res.0.8"]], "data/object/components/cluster/rna_snn_res.0.8.txt", row.names = TRUE, sep = "\t")
# # Variable features
# write.csv(VariableFeatures(so_qc), row.names=FALSE, "data/object/components/variable_features/seurat_variable_features.csv")
# # Count matrix
# writeMM(GetAssayData(so_qc, slot = "counts") %>% t(), "data/object/components/slots/seurat_counts.mtx")
# write.table(GetAssayData(so_qc, slot = "counts") %>% rownames(), row.names=FALSE, "data/object/components/slots/seurat_counts_genes.csv")
# write.table(GetAssayData(so_qc, slot = "counts") %>% colnames(), row.names=FALSE, "data/object/components/slots/seurat_counts_cellid.csv")
# # Normalized matrix
# writeMM(GetAssayData(so_qc, slot = "data") %>% t(), "data/object/components/slots/seurat_data.mtx")
# write.csv(GetAssayData(so_qc, slot = "data") %>% rownames(), row.names=FALSE, "data/object/components/slots/seurat_data_genes.csv")
# write.csv(GetAssayData(so_qc, slot = "data") %>% colnames(), row.names=FALSE, "data/object/components/slots/seurat_data_cellid.csv")
# # Scaled matrix
# writeMM(GetAssayData(so_qc, slot = "scale.data") %>% as(., "dgCMatrix") %>% t(), "data/object/components/slots/seurat_scale_data.mtx")
# write.csv(GetAssayData(so_qc, slot = "scale.data") %>% rownames(), row.names=FALSE, "data/object/components/slots/seurat_scale_data_genes.csv")
# write.csv(GetAssayData(so_qc, slot = "scale.data") %>% colnames(), row.names=FALSE, "data/object/components/slots/seurat_scale_data_cellid.csv")
# # Reductions
# write.csv(so_qc@reductions$rna_umap_nno[[]], "data/object/components/reductions/seurat_rna_umap_nno.csv")
sessionInfo()
R version 4.1.0 (2021-05-18) Platform: x86_64-conda-linux-gnu (64-bit) Running under: CentOS Linux 7 (Core) Matrix products: default BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets [8] methods base other attached packages: [1] stringr_1.4.0 patchwork_1.1.1 [3] RColorBrewer_1.1-2 ggplot2_3.3.3 [5] Matrix_1.3-4 dplyr_1.0.6 [7] celldex_1.2.0 SingleR_1.6.1 [9] SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0 [11] Biobase_2.52.0 GenomicRanges_1.44.0 [13] GenomeInfoDb_1.28.0 IRanges_2.26.0 [15] S4Vectors_0.30.0 BiocGenerics_0.38.0 [17] MatrixGenerics_1.4.0 matrixStats_0.59.0 [19] SeuratObject_4.0.1 Seurat_4.0.2 loaded via a namespace (and not attached): [1] uuid_0.1-4 AnnotationHub_3.0.0 [3] BiocFileCache_2.0.0 plyr_1.8.6 [5] igraph_1.2.6 repr_1.1.3 [7] lazyeval_0.2.2 splines_4.1.0 [9] BiocParallel_1.26.0 listenv_0.8.0 [11] scattermore_0.7 digest_0.6.27 [13] htmltools_0.5.1.1 fansi_0.5.0 [15] memoise_2.0.0 magrittr_2.0.1 [17] ScaledMatrix_1.0.0 tensor_1.5 [19] cluster_2.1.2 ROCR_1.0-11 [21] Biostrings_2.60.0 globals_0.14.0 [23] spatstat.sparse_2.0-0 colorspace_2.0-1 [25] rappdirs_0.3.3 blob_1.2.1 [27] ggrepel_0.9.1 crayon_1.4.1 [29] RCurl_1.98-1.3 jsonlite_1.7.2 [31] spatstat.data_2.1-0 survival_3.2-11 [33] zoo_1.8-9 glue_1.4.2 [35] polyclip_1.10-0 gtable_0.3.0 [37] zlibbioc_1.38.0 XVector_0.32.0 [39] leiden_0.3.8 DelayedArray_0.18.0 [41] BiocSingular_1.8.0 future.apply_1.7.0 [43] abind_1.4-5 scales_1.1.1 [45] DBI_1.1.1 miniUI_0.1.1.1 [47] Rcpp_1.0.6 viridisLite_0.4.0 [49] xtable_1.8-4 reticulate_1.20 [51] spatstat.core_2.1-2 bit_4.0.4 [53] rsvd_1.0.5 htmlwidgets_1.5.3 [55] httr_1.4.2 ellipsis_0.3.2 [57] ica_1.0-2 farver_2.1.0 [59] pkgconfig_2.0.3 dbplyr_2.1.1 [61] uwot_0.1.10 deldir_0.2-10 [63] utf8_1.2.1 labeling_0.4.2 [65] tidyselect_1.1.1 rlang_0.4.11 [67] reshape2_1.4.4 later_1.2.0 [69] AnnotationDbi_1.54.0 BiocVersion_3.13.1 [71] cachem_1.0.5 munsell_0.5.0 [73] tools_4.1.0 ExperimentHub_2.0.0 [75] generics_0.1.0 RSQLite_2.2.5 [77] ggridges_0.5.3 evaluate_0.14 [79] fastmap_1.1.0 yaml_2.2.1 [81] goftest_1.2-2 bit64_4.0.5 [83] fitdistrplus_1.1-5 purrr_0.3.4 [85] RANN_2.6.1 KEGGREST_1.32.0 [87] pbapply_1.4-3 future_1.21.0 [89] nlme_3.1-152 sparseMatrixStats_1.4.0 [91] mime_0.10 compiler_4.1.0 [93] interactiveDisplayBase_1.30.0 filelock_1.0.2 [95] curl_4.3.1 plotly_4.9.3 [97] png_0.1-7 spatstat.utils_2.1-0 [99] tibble_3.1.2 stringi_1.6.2 [101] RSpectra_0.16-0 lattice_0.20-44 [103] IRdisplay_1.0 vctrs_0.3.8 [105] pillar_1.6.1 lifecycle_1.0.0 [107] BiocManager_1.30.15 spatstat.geom_2.1-0 [109] lmtest_0.9-38 RcppAnnoy_0.0.18 [111] BiocNeighbors_1.10.0 data.table_1.14.0 [113] cowplot_1.1.1 bitops_1.0-7 [115] irlba_2.3.3 httpuv_1.6.1 [117] R6_2.5.0 promises_1.2.0.1 [119] KernSmooth_2.23-20 gridExtra_2.3 [121] parallelly_1.25.0 codetools_0.2-18 [123] MASS_7.3-54 assertthat_0.2.1 [125] withr_2.4.2 sctransform_0.3.2 [127] GenomeInfoDbData_1.2.6 mgcv_1.8-36 [129] grid_4.1.0 rpart_4.1-15 [131] beachmat_2.8.0 IRkernel_1.2 [133] tidyr_1.1.3 DelayedMatrixStats_1.14.0 [135] Cairo_1.5-12.2 Rtsne_0.15 [137] pbdZMQ_0.3-5 shiny_1.6.0 [139] base64enc_0.1-3